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ABSTRACT 

We include feedback in global hydrodynamic simulations in order to study the 
star formation properties, and gas structure and dynamics, in models of galactic 
disks. In previous work, we studied the growth of clouds and spiral substruc- 
ture due to gravitational instability. We extend these models by implementing 
feedback in gravitationally bound clouds: momentum (due to massive stars) is 
injected at a rate proportional to the star formation rate. This mechanical energy 
disperses cloud gas back into the surrounding ISM, truncating star formation in 
a given cloud, and raising the overall level of ambient turbulence. Propagating 
star formation can however occur as expanding shells collide, enhancing the den- 
sity and triggering new cloud and star formation. By controlling the momentum 
injection per massive star and the specific star formation rate in dense gas, we 
find that the negative effects of high turbulence outweigh the positive ones, and 
in net feedback reduces the fraction of dense gas and thus the overall star for- 
mation rate. The properties of the large clouds that form are not, however, very 
sensitive to feedback, with cutoff masses of a few million M©, similar to obser- 
vations. We find a relationship between the star formation rate surface density 
and the gas surface density with a power law index ~2 for our models with the 
largest dynamic range, consistent with theoretical expectations for our model of 
disk flaring. We point out that the value of the "Kennicutt-Schmidt" index found 
in numerical simulations (and likely in nature) depends on the thickness of the 
disk, and therefore a self-consistent determination must include turbulence and 
resolve the vertical structure. With our simple feedback prescription (a single 
combined star formation event per cloud), we find that global spiral patterns are 
not sustained; less correlated feedback and smaller scale turbulence appear to be 
necessary for spiral patterns to persist. 

Subject headings: galaxies: ISM - ISM: kinematics and dynamics - ISM: star 
formation - turbulence 
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Introduction 



A crucial intermediary for the formation of stars in the ISM is the gaseous cloud. Stars 
form deep within Giant Molecular Clouds (GMCs), and GMCs themselves may be embed- 
ded in larger molecular and atomic s tructures, which are referred to as giant rnolecu lar as- 
sociations (GMAs) and superclouds (IVogel et al.lll988l : lElmegreen fc Elmegreenlll983l ). The 
dispersal of cloud gas, resulting from the ionizing radiation from newly born stars, stellar 
winds, and supernovae, limits the lifetimes of GMCs and therefore determines their net star 
formation efficiencies. Supernovae (SN) also play a significant role in r naintaining and /or 
determining the thermal phase balance of the ISM (ICox fc Smith! Il974j : iMcKee fc Ostriker 
1977 : Norman fc Ikeuchi 1989 : de Avillez fc Breitschwerdtl 2004 ). and simple estima tes sug- 



gest t hat SN may be the main s ource of turbulence, at least in the diffuse ISM (e.g. ISpitzer 



19781 : iMac Low fc KlessenI |200J) . Turbulence in both the diffuse and d ense ISM is in turn 



considered one of the primary raechanisms regulating star form ation (e.g. lElmegreen fc Scalo 
2004J : iBallesteros-Paredes et al.ll2007l : iMcKee &: Ostrikerll2007l ) . Since feedback from star for- 
mation is linked to the formation, evolution, and destruction of GMCs, the overall process 
may be self-regulating. 

The formation and growth of clouds depends on the gravitational stability of the diffuse 
gaseous environment. In disk galaxies, galactic rotation and thermal pressure, among other 
factors, act to oppose the growth of self-gravitating perturbations. The Toomre Q parameter 
indicates the susceptibility of axisymmetric perturbations to grow in uniform thin disks: 
for Q < 1, the surface density is sufficiently large for gas self-gravity to overwhelm the 
restoring effects of Coriolis forces and pressure (jToomrdll964l ). Non-linear simulations have 
shown that for non-axisymmetric perturbatio ns, and including the effects of disk thickness 
and stel l ar grav ity, the threshold is Q 1.5 flKim et ahlbooa Kim k Ostriker! boOll . 120071 : 
Li et al.l l2005bl ). The observed drop off in star formation activity traced by Ha at large 
radii supports the idea that stars preferentially fo r m in gravitationally unstab le regions with 
densities above a critical value ( Kennicutt 1989 : Martin fc Kennicutt 200ll )l^ In general, 
magnetic fields cannot prevent but only slow the collapse of gas. In conjunction with other 
physical mechanisms, magnetic fields may in fact enha nce instability, as is the case when the 
magneto- rotatior ial instability (M RI, iKim et al.ll2003l ) is present, or via the magneto- Jeans 
instability fMJI. iKim et al.ll2002h . 



^Deep Ha (jFerguson et allllQQSh and UV observations |Thilker et al.ll2007l: [Boissier et al-lboOTh indicate 
that a fraction of spiral galaxies have extended outer-disk star formation, but at much lower levels than 
in inner disks. In high-Q environments, clouds may grow to become self-gravitating by successive inelastic 
collisions. 
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Star formation must commence soon after the gas acc umulates to form massive clouds, 
because almost all GMCs contain stars (IBlitz et al.l 120071 ). Ionizing radiation from newly 
formed stars subsequently dissociates the molecules, and H II region expansion disperses the 
surrounding gas; some fraction of the gas may remain molecular, but in unbound clouds. 
The massive O and B stars reach the end of their lifetimes in ~2 - 20 Myr, with those over 
8 M0 ending as SN. The cumulative effect of feedback from all the co ntiguously formiri g 
stars contributes to the short estimated GMC lifetimes of ~20 Myr (e.g. iBlitz et al.l 120071 ). 
The feedback from star formation could potentially prevent the formation of stars in nearby 
regions by driving turbulence and dispersing gas, but it could potentially trigger collapse 
events as well. Collisions between SN blast waves can result in sufficiently large densities 
for gas to collapse and form stars. It is not yet understood whether (or when) "positive" or 
"negative" feedback effects dominate; exploring this issue is one of the goals of the present 
work. 

Despite the host of processes that impact the formation of stars, observations have 
shown a clear correlation between the star formation rate density T^sfr and the gas surface 
density S, with power-law forms (in actively star- forming regions) 

ESFR OC S^+P, (1) 



now known as Kennicutt-Schmidt laws (jSchmidtlll959l : lKennicuttlll998bl ). Power law indices 
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relations have been identified for a wide range of disk galaxies at low and high redshifts. 
Both global and local versions of the T^sfr relations have been explored. In the for- 
mer, surface densities are globally averaged within some outer radius; in the latter, averages 
are over radial annuli or smaller regions. A second empirical law obtained by Kennicutt is 
TiSFR ~ O.lS/torfe, where torb is the local orbital time of the gas. 

Many theoretical studies have attempted to explain the observed relations between the 
star formation rate and the gas surface density. Simple analytic prescriptions can be obtained 
that depend on the star formation efficiency per cloud free -fall time or cloud lifetim e, and 
yield consistency with the "orbital time" empirical relations (IMcKee &: Ostrikerll2007l ). Using 
global 3D numerical simulations including g as self-gravity, a prescri ption for star formation. 



Li et al. 



Tasker fc BrvanI (120061) found power law slopes 



and feedback in the form of thermal energy, 
in T,sFR OC E^+P similar to observed values, 
found both slopes {p ~ 0.6) and normalization similar to those in iKennicuttI (Il998bl ) {p 



fl2005al . I2OO6I) using SPH simula tions 



0.4). Their simulations included gravity and sink partic l es to track the collapsing gas, but 
did not treat feedback. Recently, [Robertson fc Kravtsovl (120081 ) performed simulations that 
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included molecular cooling, and found that the power law indices obtained by fitting T^sfr oc 
j^i+P generally steeper if all of the gas, rather than just molecular gas, is included; this 
is consistent with recent observational results. 

In this work we investigate how SN driven feedback affects subsequent star formation 
in gas disks, including star formation rates. We model feedback with a direct momentum 
input, rather than using a thermal energy input (when underresolved, the latter approach 
suffers from overcooling and the resulting momentum input is too low). Our work also differs 
from other recent simulations in our approach to treating disk thickness effects; these can be 
very important to determining the star formation rate, but direct resolution requiring zones 
< 5 pc in size can be prohibitively expensive to implement in global disk models. 

The evolut ion of large gas clouds is also relevant to studies of spiral structure. In 
previous work (IShetty &: Ostrikerl |2006| . hereafter Paper I) we simulated global disks with 
an external spiral potential, and found that gravitational instability causes gas in the spiral 
arms to collapse to form clouds with masses ~ 10^ M©, similar to masses of GMAs and 
HI superclouds. We found that gas self-gravity is also crucial for the growth of spurs (or 
feathers), which are i nterarm features that are connected to the spiral arm clouds (see also 
Kim &: Ostrikerl |2002| ) . Observations have shown that spurs are indeed ubiquito us in grand 



design galaxies, and are li kely connected with large clouds in the spiral arms (lElmegreen 



1980l : iLa Vigne et al.l 120061). If grand design spiral structure i s long lasting, as hypothesized 
by density wave theory (ILin &: Shd 1 19641 : iBertin fc LinI Il996l . and references therein) , then 
feedback mechanisms dispersing the spiral arm clouds must nevertheless leave the global 
spiral pattern intact. One of the goals of this work is to assess the effect of star formation 
feedback in massive clouds on the global spiral morphology. 

Conversely, the spiral arms also affect the initial formation of clouds, therefore also 
impacting the star formation process. Observations show that most Ha emission in grand 
design galaxies occurs downstream from the primary dust lanes. An explanation for these 
observations is that gas is compressed as it flows through the spiral potential minimum, 
leading to cloud formation; then at some later time stars form within these compressed gas 
clouds. Consensus on the exact nature of spiral arm offsets has not yet been reached, however, 
owing to both observational limitations and diverse theoretical views on the star formation 
process. Fu rther, the re l ative importance of spiral arm triggering is still not completely 
understood. IVogel et al.l (119881 ) found that the star formation efficiency (in molecular gas) 
in the spiral arms of the grand design galaxy M51 is only larger by a factor of a few compared 
with interarm regions. Other observational studies comparing star formation rates in grand 
design spiral galaxi e s and those with out strong spiral structure found similar results (see 
Knapen et al.l Il996l ; iKennicuttI Il998al . and references therein). As a result, density waves 
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may primarily gather gas in the spiral arms (enhancing the ability to form GMCs), but 
may not significantly affect the star formation efficiency within any given molecular parcel. 
Without a large-scale density wave, a similar fraction of gas might still collapse (per galactic 
orbit) to form clouds via other mechanisms (including large-scale gravitational instabilities), 
but not in a coherent fashion. Here, we explore the differences in cloud formation properties 
in gaseous disks with and without an external spiral driving mechanism. 

In this work we are interested in the effect of feedback from star formation in large clouds, 
such as GMCs and GMAs, on the star formation rate, as well as on the overall dynamics 
and subsequent cloud formation in galaxy disks. This work extends the models presented in 
Paper I: numerical hydrodynamic simulations of global disks with gas self-gravity. With the 
resolutions of our models, massive GMAs do not fragment into smaller GMCs, so significant 
energy input is required to unbind the gas in these concentrations. If this energy is provided 
by star formation feedback, multiple massive stars would be needed to destroy the GMAs. 
In this work, we model feedback by considering the impact on large clouds of single energetic 
events. In practice, this could represent multiple correlated SN; this can also be considered 
simply as an expedient but cleanly parameterizable feedback model at one extreme of the 
range of event correlatioii^. We then study the resulting nature of the turbulent gaseous disk, 
as well as the formation and evolution of the clouds that form in the turbulent medium. In 
the next section, we describe our numerical simulation approach, including model parameters 
and the feedback algorithm. We then present and analyze our simulation results in ^ In 
^ we discuss our results in the context of other work, and summarize our conclusions. 



2. Modeling Method 

2.1. Basic Hydrodynamic Equations 

To study the growth and destruction of clouds in a gaseous disk, we simulate the evolu- 
tion of the gaseous component by integrating the equations of hydrodynamics. As in Paper 
I, we include the gravitational potential of the gas. Our models are two-dimensional, except 
that vertical structure of the disk is included in the calculation of self-gravity, embodied 
in a function f(z) (see below and the Appendix). The governing hydrodynamic equations, 
including self-gravity, are: 

f)y 

— + V-(Sv) = 0, (2) 



^In future work, we intend to explore how the degree of feedback correlation affects the results. 
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+ V • vv + -vn = -v($e.i + ^) - (3) 

V^^ = 4nGf{z)E. (4) 

Here, S, v, and 11 are the gas surface density, vertically averaged velocity, and vertically 
integrated pressure, respectively, and Vc is the unperturbed circular orbital velocity. For 
simplicity, we assume an isothermal equation of state, so that 11 = c^E, where Cg is the 
sound speed. The term $ is the gaseous self-gravitational potential. To grow gaseous spiral 
arms, we include an external spiral potential ^ext to model the perturbation produced by 
the non-axisymmetric stellar distribution, which is specified at time t in the inertial frame, 

by 

^ext{R, 0; t) = ^extfi cos[m0 - 0o(-R) - mflpt] (5) 

where m, (f)o{R), and ftp are the number of arms, reference phase angle, and spiral pattern 
speed, respectively. We only consider models with a constant pitch angle i, so that 

00 (^) = ln(i?) + constant. (6) 

tan i 



2.2. Model Pcirameters 

Similar to Paper I, the sound speed and rotational velocity Vc are constant in space 
and time, = 7 km s~^, and = 210 km s~^. We adopt the code unit of length Lq = 1 
kpc. Using Cs as the code unit for velocity, the time unit to = Lq/cs = 1.4 x 10® years, which 
corresponds to one orbit torb = ^tc/Qq at a fiducial radius Rq — LoVc/2nCs = 4.77 kpc. Our 
results will scale to other values of Ro and Lq with the same ratio, as well as to models with 
the same ratio Vc/cs — 30. 

In Paper I, we explored different external spiral potential strengths, 

vl tan i 

which is the ratio of the maximum radial perturbation force to the radial force responsible 
for a constant rotational velocity Vc- Wc found that spurs form in disks with strong external 
potential strengths. Since one of our objectives is to assess the evolution of the spurs in disks 
including feedback, here we only simulate disks with F — 10%, for both 2 arm and 4 arm 
spiral galaxies {m—2 and m—A). The corotation radius of 25 Lq corresponds to 25 kpc and a 
pattern speed of 8.4 km s"^ kpc~^, for spiral models using our fiducial parameters. We also 
simulate disks with no external spiral forcing. 
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In our computation of gas self-gravity, we include the effect of the thickness of the disk 
via f{z), which also acts as softening. We assume a Gaussian vertical gas distribution, 
with scale height H cc R, so the disk flares at larger radii (see Appendix). For a given 
surface dens ity, the effective niidpla ne density is given by po = S/(ifv^27r). As described in 
Paper I and iKim fc Ostrikerl (120071 ). including the effect of thickness provides an important 
stabilizing effect on the disk. For most of our simulations, we use H/R = 0.01. 

As in Paper I, the Toomre parameter Qo = kqCs/ItcGTio) and the surface density So at 
Rq are related by: 



GLqQo Qo 



32 

— Mopc- 



7kms 



kpc 



For fiat rotation curves, the epicyclic frequency k, = ^/2^l = \^Vc/ R- Our models initially 
have S oc R'^, so that Q is constant for the whole disk. 



2.3. Numerical Methods 



Since this work is an extension of previous work, we refer t he reader to Paper I for a 
description of the cylindrical-symmetry version of the ZEUS co de Jstone fc Normanlll992al jbl) 
that we use to carry out our simulations. We use a parallelized version of the hydrodynamic 
code and gravitational potential solver, allowing us to increase the number of zones in the 
grid relative to the models of Paper I. For our standard grid we set the azimuthal range 
to - 71/2 radians and the radial range to 4 - 11 kpc. We implement outflow and periodic 
boundary conditions in the radial and azimuthal directions, respectively. These models have 
1024 radial and 1024 azimuthal zones. Since the radial grid spacing is logarithmic, the 
resolution varies: the linear resolution in each dimension {AR, i?A0) varies from ~(4 pc, 6 
pc) in the innermost region t o ^(11 pc, 17 pc) at t he outer boundary. These high resolutions 
allow the Truelove criterion (ITruelove et al.l 119971 ) to be satisfied throughout the simulation 
as gas collapses to form self gravitating clumps. 

In this work we use a different method to compute the gravitatio nal potential fr om 
that in Paper I. Here, we use a method derived fro m that described by iKalnajd (jl97ll ) in 
polar coordinates (see also lBinney fc Tremaindll987l ). This method employs the convolution 
theorem for a disk decomposed into logarithmic spiral arcs. We implement softening to 
account for the non-zero thickness of the disk. We describe the method in detail in the 
Appendix. 



We note that for simulations with the standard grid and including a spiral potential, 
the limit in azimuth requires that m=4 (4 arms). However, we also explore some models 
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with m=2 patterns, with the azimuthal range - vr, using twice as many azimuthal zones 
than the standard grid so that the physical resolution of both simulations are equivalent. 



2.4. Feedback: Event Description and Algorithm 

Equations ©-dlD only describe the flow as gas responds to self gravity, and to the 
external spiral perturbation, if one is present. However, those equations do not describe any 
feedback that would occur after a self-gravitating cloud forms and fragments into smaller- 
scale structures, ultimately forming stars with a range of masses. In the real ISM of galaxies, 
clouds are dispersed by the combination of photo-evaporation by UV radiation from massive 
stars, and the "mechanical" destruction by expanding HII regions and SN. 

We include in our simulations a very simple feedback prescription by implementing 
"feedback events," each representing momentum input from a number of SN (or, alterna- 
tively, multiple overlapping expanding HII regions). The specific SN rate, Rsn, averaged 
over all mass Mdense above a chosen threshold density in a galaxy is 

Number of Supernovae , , 

Mdense " time 

where Nsn is the number of supernovae. When this rate is applied to an individual dense 
cloud of mass M^i with a lifetime td, the average number of SN in the cloud will be 

Nsn = Rsn ■ ■ U. (10) 

If the total mass of stars of all masses formed per single SN is Msn, and the star formation 
efficiency over a cloud lifetime is esF, then 

Nsn = ^SF^. (11) 

Msn 

Equating expressions ffTOj) and ffTTj) . the mean cloud lifetime is 

^SNJ-^J-SN 

In a given time interval 5t, such as the time between successive computations in the 
numerical evolution, the probability P that a cloud (of mean lifetime td) is destroyed is 
St/tci- Thus, 

P = 5t- RsNMsN/esF (13) 

In our algorithm, clouds are defined as regions above a chosen density threshold. If a par- 
ticular zone is a local density maximum, that zone is selected as the center of the feedback 
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event. For any such identified cloud, a star formation event is initiated witli a probability 
per timestep given by equation f|T3|) . In each cloud that is determined to undergo feedback, 
gas is evenly spread out in a circular region with a prescribed bounding radius. Gas in 
each zone in the circular region is assigned an outward velocity (relative to the center) to 
expand the feedback "bubble." A constant azimuthal velocity is also added such that total 
galactocentric angular momentum is conserved. The velocity profile inside the bubble is 
proportional to the distance from the bubble center. For most of our simulations, we choose 
the radius of the feedback bubble to be 75 pc, which corresponds to 12-23 pixels, depending 
on the radial location. In this way, the initially collapsing cloud gas is forced back into the 
surrounding ISM. 

In our simulation, we only consider the isothermal expansion of the clouds, since we 
assume an isothermal equation of state. Thus, we can only consider the net energy input at 
a stage when expansion of the shell has become strongly radiative. Numerical simulations 
show that for a single SN of energy Esn ~ 10 ^^ ergs, the radial momentum d uring the 



radiative stage is Prad ~ 3 - 5 x 10^ Mq km s~ ^ (I Chevahei] 1 1 9 74 ICioffi et al.lll988l ). During 
the subsequent evolution of the bubble, the shell momentum Pgh is conserved, and is equal 
to Prad- Wind-driven and pressure-driven H II region bubbles similarly are accelerated to 
reach a final momentum Prad- 

For a total number of massive stars formed given by equation ffTTj) . and assuming cor- 
relation in time, the total momentum applied to the shell is 

Psh = NsNPrad = ^SF-jTJ Prad (14) 

MsN 

The shell velocity Vsh is 

= Psh/M,, = ^SF§^§^- (15) 

Here, Mg^ is the sum of and any ambient gas in the (circular) feedback region. Assuming 
MsN = 100 Mq and egp = 0.05, for Prad = 3 x 10^ Mq km s-\ Vsh = 150 km s-^x (M^^/M,,,). 

Given our feedback prescription, the two key parameters are the probability per unit 
time for cloud destruction (eq. [I3]), and the momentum input in the feedback event (eq. 
[13]). For the simulations presented here, we explore a range in the rate Rsn and in the 
momentum input per massive star, Prad- The specific SN rate is set either to Rsn = (10^ 
Mq X 50 yr)"^ = 2 x 10~^^ Mq'^jt'^ (comparable to that in the Milky Way), or ten 
times that rate (these models are denoted by Rsn = 1 or 10 in Table [1]). Since Msn/^sf 
in equation (ITHl) appears as its inverse in equation (HM. we fix M sn = 100 Mq for all 



simulations, motivated by the initial mass function of iKroupal (120011 ) . and explore variations 
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in esF- Scaling to fiducial values, we then have 

(16) 



2 X 10-11 m;;^ 



yr" 



for the typical lifetime of dense clouds. The momentum Prad is set either to 3.4x10^ or 
3.4x10^ Mq km s~^, in order to allow for a range in feedback energy and out-of-plane losses 
(venting from the galaxy) that reduce Prad for a given energy input. 

We note that for low values of Prad, the energy input will not be sufficient to destroy a 
dense, bound cloud. In particular, a cloud of surface density T^d will become unbound only 
if Vsh ^ (2G)i/2(7rSciMd)^/^ For = 200 Mopc'^ and M^i = 10^ M©, the minimum shell 
velocity is ~ 15km s~^. For our larger value of Prad, this inequality is comfortably satisfied, 
but for the smaller value it is not. We indeed find that for the low Prad models, clouds are 
not destroyed by feedback. For these models, then, the ratio td = esp/iRsNMsN) becomes 
the mean interval between (non-destructive) feedback events in a given cloud.. 

Before any feedback, the spiral models are executed for some time to allow gas to 
concentrate (due to self gravity) and form clouds in spiral arms. In simulations without 
spiral forcing, condensations begin to grow due to an initial 0.1% density perturbation. As 
a result of shear, the first structures that form are large scale fiocculent spiral-like features, 
which we termed "sheared background features" in Paper I. Gas in these features then 
collapses to form distinct clouds. Thus, we wait until some threshold density is reached 
before feedback occurs. For most models, the threshold density is S/Sq = 10 (this sets the 
threshold S at 320 Mq pc'^ for Qo=l, and 160 Mq pc'^ for Qo=2). We hereafter refer to 
any contiguous structure in our simulation with a density above this chosen threshold as a 
"cloud," regardless of whether the given structure hosts a feedback event or not. 

We note that with our feedback prescription, the star formation rate is given by the 
mass in dense gas (i.e. exceeding the threshold surface density) times RsnMsn- This 
linear relation is supported by the shallo w slopes of S.sfr versus T^mni (as observed i n CO 



emission) . In some other rece nt work (e.g. iLi et al.ll2005al . l2006l : iTasker fc Bryanll2006l . 12008 



Robertson fc Kravtsov 20081 ). the star formation rate is taken as equal to the mass (with 



density above some threshold) divided by the free-fall time at that density, times some 
efficiency factor. Our prescription is therefore equivalent to choosing a ratio of efficiency 
over free-fall time at the surface density threshold of e///^// = RsnMsn = (5 x 10^^^yr)^i. 
Since the mean internal density within real GMCs (which have surface densities similar to 
our critical threshold) is ~ 100 cm"'^, with corresponding free-fall time of 4 Myr, our models 
would cover a range of star formation efficiencies per free-fall time of e/j ~ 1 — 10%. 
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3. Simulation Results 

We first present simulations witli standard grid parameters, without spiral structure. We 
then show results of simulations including spiral structure, as well as simulations pertaining 
to different radial regions. 

Table [1] shows the initial conditions of the standard set of models we present, as well 
as the relevant parameters controlling the feedback events. Column (1) lists each model. 
Column (2) shows the initial Toomre Q parameter which is initially constant for the whole 
disk. Column (3) indicates the number of arms, all with F=10%. Column (4) gives the SN 
rate, which is required for setting the probability that a feedback event occurs in a cloud 
(see eq. [I3]). Column (5) shows the assumed star formation efficiency, and column (6) gives 
the adopted momentum input per massive star. For these models H/R = 0.01. 

3.1. Disks without Spiral Structure 

Figure [T] shows a snapshot of model with Qo=i, at time t/torb = 0.84, without an 
external spiral potential and before any feedback. As discussed in §2.41 trailing features 
grow due to the self-gravity and shear in the disk (see Paper I for details). The most dense 
structures grow as sheared, trailing features. It is in these regions where the first SN will 
occur to disperse the dense gas. 

Figure [2] shows a snapshot of model QIA, at time t/torb = 1.125. For model QlA, 
the SN parameters are all at the low end of the range. At the time of this snapshot, 105 
feedback events have occurred, in clouds which have mean Md = 1.2 x 10^ M0. The main 
difference between Figures [2] and [T] is the shape of the trailing features. The feedback events 
have caused the features to become fragmented at some locations. However, dispersal of 
gas due to feedback was not sufficient to prevent or reverse the inflow of gas into the high 
density agglomerations. Either the SN do not occur rapidly enough, or do not have enough 
momentum to alter the basic morphology. Even increasing both the SN rate by a factor of 10, 
and doubling the star formation efficiency makes little difference; the strong self-gravitational 
force from the trailing features keeps much of the gas in those structures. Increasing the SN 
momentum (or equivalently the velocity) by up to a factor of 8 still does not significantly 
affect the outcome: much of the gas is contained in the sheared structures at any given time. 

It is only when P^ad is increased to 3.4x10^ Mq km s~^, along with increasing Rsn by a 
factor of 10 and esF to 0.05, that we find a significant difference compared to the case QlA, 
as in Model QID shown in Figure [31 The velocity is sufficiently large to drive gas away from 
the density maxima of the trailing structures. Further, the rate is high enough that a large 
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number of events occur to significantly alter the morphology, in comparison with Figure [2j 
Feedback events in this model are so frequent and energetic that collisions between bubbles 
occur. In some instances, such collisions create density enhancements that later result in 
more collapse and subsequent feedback along the bubble interface. At time t/torh = 1-125 
(Fig. [3](a)), we can still make out the underlying loci of the initial structures formed by 
gravitational instability and shear, though 537 feedback events have occurred up to this 
time. Yet, after an additional 26 Myr and 75 feedback events (Fig. [3](b)), the dominant 
large scale features do not have a single pitch angle. Further, the locations of many of the 
bubbles are clustered. Though gas is driven away from the initial structures formed before 
feedback, at later times clouds form in clusters near the initial density maxima, and where 
feedback bubbles overlap. Qualitatively, the features in the di sk, consisting of filament s and 



bubbles, are similar to the global models including feedback of IWada fc NormanI (|2001| ). We 
discuss the masses of the clouds in both "non-spiral" and spiral models in the next section 

(M). 

In disks with Qo = 2, sheared features will also grow, but need more time to develop 
than in the Qq = 1 disks. Due to its relative stability, after t/torb = 2 only a few clouds have 
formed. As a result, implementing feedback does not affect the majority of the disk. To 
study the effect of feedback in Qo = 2 disks, another mechanism is necessary to grow clouds 
everywhere in the disk. We thus simulate Qo=2 disks with an external spiral potential, and 
then implement feedback to destroy the spiral arm clouds that grow. 



3.2. Disks with Spiral Structure 

In disks with spiral structure, the stellar spiral potential acts as a source of perturbation; 
the compression of gas as it flows through the potential eventually leads to the growth of 
self gravitating clouds. We explore the effect of feedback on the morphology of the gaseous 
spiral arms and interarm spurs, as well as any subsequent cloud formation. 

Figures miS] show snapshots of models with m = 4, for Qq = 1 and Qq = 2, without any 
feedback. In the spiral models, the growth of spiral arm clouds occur sooner than clouds 
formed by natural instabilities in a rotating self-gravitating diskjfl Figures HJlS] show snapshots 
of models without any feedback, though in the Q = 1 snapshot (Fig. Hj) the densities have 
surpassed the threshold density S/Sq = 10 chosen for models with feedback. Note that while 



•^In our models, the amplitude of the spiral perturbation is "turned-on" gradually, reaching the maximum 
amplitude F at t/torb = 1- Due to this imposed "turn-on" time, the growth rate of GMCs in our models is 
not representative of actual GMC formation timescales. 
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the Q = 1 model (with strong self-gravity) shows dense gas knots within the arm, the Q = 2 
model (with weaker self-gravity) shows spur-like features; gas does not collapse as promptly. 

For both models QISA and Q2SA, the SN momenta are insufficient to offset the growth 
of clouds and spurs resulting from the spiral potential. After a feedback event, the dispersed 
cloud gas flows back toward the spiral arm. As a result, clouds continue to grow over 
time. Further, the spurs also continue to grow in density. Without feedback, self-gravitating 
spiral arm clouds cause the surrounding gas to flow in with large velocities. Eventually, 
the simulations have to be stopped because the Courant time is too small. The time when 
the simulation ceases, depending partly on our choice of the minimum acceptable Courant 
time, also depends on which clouds are (randomly) selected for feedback; clouds that have 
produced large inflow velocities would have to be dispersed for the simulation to continue to 
evolve. 

We again find that large SN momenta are required to sufficiently disperse clouds so that 
immediate re-collapse does not occur. For such models, the SN rate has an effect on the 
number of subsequent clumps formed. FigureElshows a snapshot of models QISC and QISD, 
~ 21 Myr after the first feedback events. At this time, 53 feedback events have occurred 
in model QISC, and 540 in model QISD. In model QISC, it is clear that most, if not all, 
feedback events originated in the spiral arms. However, in model QISD, many feedback 
events have occurred in interarm regions. The spiral arms are not as identifiable, though at 
this time the remnants of spurs are still identifiable. Further, model QISD contains many 
more clumps than model QISC. The enhanced SN rate has caused the collision of more shell 
remnants, which lead to formation of self gravitating clumps at the interfaces. In both cases, 
feedback events have caused gas to be dispersed from the arms, eventually removing any 
trace of the underlying spiral potential, as can be seen in Figure [71 

Figure [S] shows the histogram of the masses of the clouds that hosted feedback 
eventqj in models QlA, QID, QISA and QISD. In all cases, the maximum mass of the 
clouds is below 10'' Mq, and the means and medians for the distributions lie in the range 
0.5 — 2.2 X 10^ Mq. In model QlA, most feedback events have occurred in the large scale 
sheared features that grow due to gravitational instability. However, in model QID, some 
fraction of the feedback events have occurred in regions of colliding fiows. The histogram 
suggests that clouds formed by colliding fiows have characteristically lower masses than those 
formed in the large scale sheared features. Similarly, in model QISA, most feedback events 
have occurred in the spiral arms, since most clouds form in the arms. On average, the clouds 



*From equation feedback events in model QlA and QID on average consist of 300 and 350 SN, 

respectively. 
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in model QISD have lower masses, with many clouds formed due to colliding flows initiated 
in earlier feedback events. For a power law in the mass distribution, dN/dlogM oc M~", 
the distribution in the high end masses for model QISD (below the cutoff at log(M) = 6.4) 
gives a ~ 0.6. This slope and the upper limit in cloud n aasses is similar to the r ange 
and the upper limit in the observed masses of GMCs (see iMcKee fc Ostrikerl 120071 . and 
references therein). A histogram of the masses of all clouds at any given time, normalized 
by the correct probability St/td, reproduces the overall shape of the histogram of the of 
masses of clouds with feedback. A detailed analysis of the cloud mass distribution is not 
appropriate here because many of the lower mass clouds are not well resolved, and because 
higher-mass clouds would be subject to turbulent fragmentation that we cannot follow. 
Higher resolution simulations are therefore required to obtain more complete cloud mass 
distributions. Nevertheless, it is clear that the upper mass limits for clouds in all models are 
similar to those in real spiral galaxies. 



3.3. Star Formation Properties 

3.3.1. Star Formation Rates and Turbulence 

For comparison to observations, two quantities of interest are the star formation rate 
SFR and the turbulent velocity Vturb- In each simulation, we record each feedback event to 
determine the SFR. For some chosen time bin At, we compute 

SFR = esF^^, (17) 

where ^ M^i is the total mass of all gas in clumps (i.e. above the chosen threshold surface 
density) that have undergone feedback events in the chosen time interval. (Recall that the 
mean lifetime of clouds, or the mean interval between star formation events if they are 
non-destructive, is given by eq. [12].) 



We define the turbulent velocity as the RMS sum of any non-circular velocities, weighted 
by the corresponding mass: 



'^turb 



where Ai^j is the area of each zone and only non-circular velocity components are considered: 
5v = V — Vc(f). Figure [9] shows the star formation rate and turbulent velocity as a function of 
time, for the Qo=l models without spiral structure. The time bin At for our SFR calculation 
is 3 Myr. In these models, the first feedback events occur at time ~125 Myr. However, for 
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the first ~25 Myr after feedback begins, the SFR for all models is only a few M0 yr^^. 
Only ~25 Myr after the first feedback events does the SFR substantially increase, owing to 
"propagating" star formation. Further, the QID model with large feedback momenta {Prad 
= 3.4x10^ M0 km s"^) and large SN rate {Rsn = 10) has the SFR increase to ~10 M0 yr'^ 
This occurs because with large velocities and a high global rate, adjacent shells collide and 
more clouds are formed in the interfaces, which may subsequently undergo star formation. 

The bottom panel of Figure [3 shows Vturb, for all feedback models without spiral struc- 
ture, together with results from a simulation without any feedback. For the later case, we 
just allow self-gravity to grow clouds indefinitely. When we compute Vturb in the model 
without feedback considering only the low density gas, we obtain similar values. This sug- 
gests that, before any feedback, lar ge-scale motions from disk self-gravity and shear are the 



primary sources of turbulence (see iKim &: Ostrikerl 120071 ). The models with low feedback 



momentum continue the trend of Vturb established by the no-feedback case. In a few instances 
of enhanced feedback, there is a corresponding jump in Vturb- The enhanced SFR at later 
times for the model with large SN momenta also increases levels of Vturb- 

Figures [TD] and [TT] show the SFR and Vturb for the spiral models with Qq = 1 and 
Qo = 2. Comparing Figures [TO] (with spiral structure) and [9] (without spiral structure), the 
star formation rate is consistent to within a factor of 2, although slightly larger in some of the 
spiral models. The general trends from the models without spiral structure are reproduced 
in Figures [TD] and [TT] Earlier times are shown in Figure [TD], since the spiral arms cause gas 
to collapse into clouds sooner. It is clear that only in models with large SN shell velocities 
- and regardless of the input rate Rsn - do the turbulent velocities increase appreciably; 
otherwise, the turbulent velocity (as we have defined it) is dominated by effects from gas 
self-gravity. 

It is interesting to compare results from pairs of models in which one parameter is 
varied and the others are controlled. Comparing models QISE and QISB, both have the 
same esF = 0.05 and Prad = 3.4 x 10^ Mq km s~^, but the former has Rsn larger by a factor 
10. The measured SFR in QISE is a factor ~ 10 larger than that in QISB, consistent with the 
naive expectation that SFR oc Rsn- However, when we compare QISD with QISC, which 
again differ in Rsn by a factor 10, we find SFR ratios differing only by a factor ~ 4. This same 
trend is also true for models Q2SD and Q2SC. The reason for this difference in dependence 
on Rsn is that the E and B models have low Prad and low turbulence levels, whereas the C 
and D models have higher Prad and turbulence. Thus, stronger feedback causes the scaling of 
SFR to depart from SFR oc Rsn- We note that since SFR = RsNMsN^dense by definition, 
the ratios of specific SFR between any two models differ by their ratios of Rs nM dense / Mtot- 
Thus, if SFR increases at a rate less than oc Rsn, it implies that increasing Rsn decreases 
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the dense gas fraction Mdense/Mtot- 

We can directly investigate the effect of turbulence by comparing the pair QISB and 
QISC, which have the same esF = 0.05 and Rsn, but momentum input parameters differing 
by a factor 10. As noted above, this increases the turbulence level in QISC compared 
to QISB by about 10 km s~^. It also reduces the SFR in QISC compared to QISB, by 
a factor ~ 2 — 4. Similarly, QISE has lower P^ad than QISD, and a substantially lower 
turbulence level. For this pair, too, the SFR in the lower-turbulence model is higher by a 
factor ~ 3 — 5. As discussed in §11 in principle turbulence could both enhance star formation 
(by creating more dense gas in compressions), and suppress star formation (by destroying 
overdense structures with rarefactions and shear flows). Examining the evolution of QISE 
indeed shows that feedback events only slightly expand clouds, and collapse subsequently 
resumes. On the other hand, clouds in model QISD are completely destroyed after a single 
feedback event. Evidently, in the models with strong feedback-driven turbulence, the rate 
of new cloud formation from shell collisions does not compensate for the truncation of star 
formation when a given cloud is destroyed. 

The comparisons of (QISB, QISC) and (QISE, QISD) indicate that in net, the increase 
of turbulence reduces star formationjf Since the specific SFR is proportional to the dense 
gas fraction if Rsn is held fixed, these results imply that the dense gas mass fraction is lower 
when the turbulence level is higher. 

We show the relationship between the mass weighted turbulent velocity and the surface 
density in Figure [121 Most feedback events occur in high density regions. In the higher 
density regions, the difference in turbulent velocities (or the velocity dispersions) between 
models with Prad = 3.4x10^ Mq km s"^ and Prad = 3.4x10^ Mq km s"^ is > 7 km s'^. 
At lower density regions, where there have been fewer feedback events, the dispersions of all 
models are comparable. 

Figure [T3l shows the turbulent power spectrum (power oc v^) of model QID. The power is 
shown at constant wavenumbers and k^. The slopes of the power spectra range from -2.5 
to -3. For models that evolve for significant amounts of time, such as model QID, the power 
spectra are relatively independent of time. These results are consistent with turbulence 
dominated by numerous shocks, or Burgers turbulence. From Figure [131 the amplitudes of 
turbulence evidently decrease at smaller scales. 

The total turbulent amplitudes shown in Figs. lUllTT] represent the velocity dispersion 



^We note that models Q2SB and Q2SC also show the same generic behavior, but with a smaller difference 
in the turbulence level, the suppression of star formation is also lesser. 
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averaged over the whole disk. For the purposes of assessing turbulent contributions to local 
disk stability, however, only the level of turbulence within a Jeans length ~ c^/((jS) is 
relevant. Furthermore, local observations of turbulence within the Milky Way generally 
measure velocity dispersions on scales less than the disk thickness. Thus, it is useful to 
estimate the turbulent amplitudes at smaller scales than the whole disk. We do this by 
running a window (or "beam") of 1 kpc or 100 pc over the map, and finding the dispersion 
of the velocity within this window at locations separated by the window size. When all 
zones within the window are weighted equally (as is true for the velocity power spectrum), 
we find that the mean velocity dispersions for model QID on scales of 1 kpc and 100 pc are 
18 km s~^ and 6.5 km s~^, respectively. When we weight by mass, the respective velocity 
dispersions are 31 km s~^ and 10 km s~^. The larger values obtained when weighting by 
mass are indicative of the importance of dense expanding shells in driving the turbulence. 

Since turbulence adds to the total momentum flux (the ram pressure acts similarly to 
the thermal pressure), a common assumption is that the sound speed Cg can be replaced by 

clff = cl + al (19) 

in the dispersion relations that characterize stability to axisymmetric modes, where ctr is the 
radial component of the velocity dispersion. For models QISC and QISD, which have high 
Pradi we find that the mean values of aji on kpc scales are 17 and 18 km s~^, respectively. 
For the corresponding models QISB and QISE that have low Prad, on the other hand, the 
values of aR on kpc scales are 10 and 8 km s~^, respectively. Thus, the values of Cg// exceed 
Cs by a factor 1.6 for the low-turbulence models, whereas this increases to a factor 2.7 for 
the high-turbulence models. Our results discussed above indicate a decrease in the star 
formation rate with increasing Cg//; we discuss theoretical ideas related to this finding in 
§3.3.31 below. 



3.3.2. Kennicutt- Schmidt Law 

Figure [H] shows the local star formation rate per area as a function of mean surface 
density. To obtain these points, simulation data were binned in radius and time, of widths 
1 kpc and 18 Myr, respectively. Only models with a sufficient number of points, which is 
dependent to some degree on the number of feedback events, are shown. Best fit lines to the 
data points are also shown. The rates show considerable scatter, both between models with 
different parameters, as well as among points from a given model. However, where a large 
dynamic range is available, as is the case for the QID model and its extension to smaller 
radii (see below), a power law relation T^sfr S^+p is quite clear. 
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The Qo = 1 models, both with and without a spiral perturbation, and with different 
feedback parameters, generally give slopes 1+p ~ 1—3. Most of the Qq = 1 models evolve for 
sufficiently long times that gas in the first clouds that are formed are allowed to be recycled 
into subsequently formed clouds several times. The Qq = 2 models, on the other hand, give 
a variety of slopes, and the relationship between the star formation rate and surface density 
is not as well correlated as in the Qo = 1 models. For the Qo = 2 models, the number of 
feedback events is insufficient to affect much of the disk. As a result, some clouds continue 
to collapse, and the Courant condition would demand an extremely small time step; at this 
point, we halt the simulation. Since the stochastic feedback events do not result in developed 
turbulence and a steady state is not approached in the Qo = 2 models, the SFR as computed 
is sensitive to model parameters governing the feedback events. 

For some of our models, we have also run simulations of the inner regions of disks, with 
radial extent R G 0.8 — 2.2 kpc. The other parameters are the same as for the standard 
models. The only difference here, besides the radial range, is the initial surface density. 
Since Qo is constant, and Sq oc {QoR)~^, the initial surface density at all radii is increased 
by a factor 5 compared to the standard models with i? G 4 — 11 kpc. 

Figure [15] shows the star formation rate as a function of surface density for model QID 
together with the corresponding inner region model. The larger surface density does indeed 
lead to higher star formation rates, with a slope 1 + p = 2.2 that is similar to the value 
1 + p = 2.4 of the standard model. We find similar trends for other inner disk models in 
compariso n with the corresponding standard models. For comparison. Figure [TS] also shows 



data from iKennicuttI (jl998bl ). Each point indicates the globally averaged star formation rate 
for individual galaxies or their central regions (for starbursts). Though there is less scatter 
in the simulation points, the slope of the "Esfr - S relation from the simulations (~2.3) is 
larger than the slope from observational data (1+p ~1.4). At the low S end, the model 
results overlap with the observed points. 



3.3.3. Predicting Star Formation Times 

The star formation (or gas depletion) time tsp for the whole gaseous component of a 
galaxy is the time required for all the gas to be converted to stars if the star formation 
proceeds as it has been during a given interval At: 

where Mtot is the total mass in a given annulus. This quantity can be measured in our 
simulations; the summation in equation ( !20|) is taken over all clouds in which a feedback 
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event has occurred, as in equation ([T 

Observationally, if the SN rate per dense gas mass Rsn is known, the star formation (or 
gas depletion) time can also be estimated based on the total amount of gas and the portion 
in dense clouds as: 

t' (21) 

where Mdense is the total mass of gas above some chosen threshold density. Since RsnMsn = 
esp/tci from equation (|T2l) . the results of equations (|20|) and fl2Tl) should agree on average. 
With our two parameter choices Rsn = 1 or 10 (in units 2 xlO^^^ M0^^ yr^^), this implies 
t'gp = (1 or 0.1) X (Mtot/Mdense) X 5 X 10^ yr. We note that if the star formation or 
gas depletion time were computed only for dense gas (with local surface density ^ 200 
Mq pc~^), then for our prescription it would simply be equal to a constant, t'gp{dense) = 
{MsnRsn)~^ = 5 X lO'^ or 5 X 10^ yr for Rsn = 10 or 1, respectively. 

Figure [TBI shows the star formation time in different radial annuli for Model QID, as a 
function of S. The actual times, shown by the filled symbols, are computed using equation 
f l20|) . after binning the simulation data in radii of 1 kpc widths and in time with t/torb = 
0.125 widths. The open symbols show the predicted times by applying equation ( l2TI) on the 
same binned data. The predicted times agree well with the actual times. We find similar 
agreement with all other models]^ We also tested the correlation between tsF and tort, and 
found no strong correlation. This lack of correlation occurs because at later times the surface 
density profile no longer resembles the initial R^^ profile. 

What is expected, on theoretical grounds, for the value of the star formation time? 
Consider the case in which gas cycles between diffuse (gravitationally unbound) and dense 
(gravitationally bound) components. The diffuse component forms dense clouds at a rate 
Mdiff/tdiff, and the dense clouds are returned to the diffuse component plus stars over a 
cloud lifetime at a rate M dense/ id- Here Mdifj and M dense are the total diffuse and dense gas 
masses in an annulus, with corresponding surface densities when averaged over the area of 
T^diff and T^dense (the latter is not to be confused with the surface density of an individual 
dense cloud, which is much higher). Similarly, S is the surface density corresponding to the 
total mass of all the gas Mtot in an annulus. In equilibrium, the rates into and out of the 
dense component are equal, so that the star formation rate per unit area averaged over the 
annulus is 

SsFR = ^SF- — — = ^SF- = ^SF — ' = RsNMsN^dense- (22) 

idiff + tcl tdiff tcl 



®Most other models do not run for as long as the D models that have high feedback rates, because some 
dense clumps continue to collapse without feedback, eventually causing the Courant condition to be violated. 
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Including all the gas, the star formation timescale using the definition of equation fl2T]) (and 
dropping the prime) is then 

, tdiff + tcl tdiff Mtot tcl Mtot 

tsF = — = — — X — = X — . (23) 

^SF ^SF Mdiff €SF Mdense 

Since td/esF = {RsnMsn)^^ = 5 x 10^ or 5 x 10^ yr is held constant within any given model, 
the star formation time for all gas in an annulus is inversely proportional to the fraction of 
the gas above the density threshold in that annulus. If most of the gas is diffuse (as is true 
in our simulations), then tdifj ^ and tsF ~ tdiff/^SF] the star formation time is set by 
the typical time required for diffuse gas to collect into bound clouds. 

What characteristic values might be predicted for the cloud formation timescale, tdi//? 
The shortest possible timescale would be that associated with the fastest-growing Jeans 
modes in a disk. For a disk with semi-thickness H and sound speed c^, the approximate 



dispersion relation for in-plane modes is uo"^ = k'^c^ — 27rGE|/c|/(l + \k\H) (IKim et al.l 12002 
Paper I). For the fastest-growing modes (which satisfy d\uj^\/dk = 0) and for H < Cg/(7rGS) 
(i.e. less than the thickness of an isothermal disk bound only by its own gravity), the inverse 
of the growth rate is 0.3 — O.Stj, where tj = Cs/{GT,) is the thin-disk Jeans length divided 
by Cg. In reality, rotation, shear, and turbulence must all affect the cloud growth timescale 
(see below), but the Jeans time nevertheless provides a useful reference value. 

Another reference value for a structure formation timescale that is frequently used is 
the free-fall time, t// = (37r/32Gp)^. If the surface density and volume density are related 
via S = pH\/2m (as for a Gaussian density distribution), then tfj = (3v^7r^^^-f^/32G'S) 2 . 
For our "thick-disk" Poisson solver, H ex. R is adopted, so that tjf oc (-R/S)^/^. Our initial 
profiles follow R oc so that in the initial conditions t// oc S^^ oc tj. In particular, for 
the Q = 1 case, t// = 0.3tj everywhere initially. Over time, however, the surface density is 
spatially rearranged, so that the values of tj and t// are no longer strictly proportional. 

Figure [17] shows the relationships between the star formation time and the reference 
values tj and tff. While a clear correlation is evident for both relations, we find that there 
is less scatter in the tsF ~ tj relation than in tsF -~ tff relation. Further, many of the data 
points are consistent with a linear relationship tsF = 7tj, as indicated in the figure. If 
we compare to the prediction tsF = tdiff/^SF and substitute the value esF = 0.05 used in 
model QID, this yields tdiff = 0.35tj, which agrees with the simple estimate described above 
based on self-gravitating instabilities in thick disks. This result suggests that, provided the 
efficiencies of star formation in GMCs are constant and the disk is dominated by diffuse 
gas, the Jeans time in the diffuse gas controls the rate of star formation. While this result 
is quite intriguing, a high djTiamic range in a wider range of disk models is necessary to 
further investigate this relationship. 
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We note that in the dispersion relation used to predict tdijf ^ tj = Cs/ (GS), no account 
was made for turbulence. As discussed in §3.3. H the simplest phenomenological modification 
of this relation would simply be to substitute Cg Ce/f (see eq. [IH])- The results presented 
in §3.3.11 which compare SFRs for model pairs with low and high Prad, and hence different 
Ceff, are at least semi-quantitatively in support of this prescription for modifying tj. There, 
we found that an increase of Ce// by a factor of ~ 2 is associated with a decrease in the SFR 
by a factor ~3. However, the current models are not sufficient for a definitive statement. An 
important objective for future work is to test the relation between tsF and the turbulence 
level using a more extensive set of models; the velocity dispersion can be varied by tuning 
the parameter Prad- A fundamental understanding of star formation in molecular-dominated 
regions of galaxies (where the thermal velocity dispersion is dwarfed by the turbulent value) 
will depend on such investigations. 

Modeling truly three-dimensional disks, with the vertical dimension fully resolved, would 
allow for a more complete study of the correlations between tsp and the two gravitational 
times, tj and t//. Depending on the regime, vertical hydrostatic equilibrium (for an isother- 
mal medium) may be in the limit dominated by (a) the disk's gaseous self-gravity, so that 
the effective thickness of the ISM is S/(2po) = cl/^nGT,), or (b), the disk's stellar grav- 
ity, so that the effective thickness is S/(2po) = Cs<J^/ {2^/^TGT,^) oc {Q^/Q)cj./ (nGT,). Here, 
cr* and T,^: are the stellar vertical velocity dispersion and surface density, respectively, and 
is the Toomre parameter for the stellar disk. Using these two forms, if gas dominates 
the vertical gravity, then tff oc tj, whereas if the stars dominate the vertical gravity, then 
tff (X tj{Q^/Qy^'^. If galaxies evolve such that Q^^/Q is constant, then tj oc tff in either 
case; it would then be empirically difficult to establish whether tj or tff is more fundamen- 
tal for determining the star formation time. With explicit three dimensional models, on the 
other hand, it will be possible to study the dependence of tsF on tj and tff separately, with 
Qtf/Q a. tunable parameter. This represents a very interesting avenue for future research. 

4. Discussion and Summary 

4.1. Kennicutt-Schmidt Law in Simulations 

The prescription we adopt for star formation in this paper implies a constant relation 
between the mass (or mean surface density) of dense gas and the rate (or mean surface 
density) of star formation, T^sfr = RsNMsN^dense- Using this prescription, we then test 
how the star formation rate scales with the surface density of all the gas. We find that 
our simulations are consistent with scalings T^sfr oc E^"*"^ for a range of power law indices, 
but with significant scatter. In part, both the range of indices and the scatter in many 
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of our models may arise from transient effects, rather tfian describing the behavior in a 
fully-developed star-forming disk. Our simulations suggest that measured star formation 
properties are subject to transient effects, thus for meaningful theoretical predictions it is 
necessary for systems to evolve well beyond the initial state. 

For our strong-feedback model that most closely reaches an equilibrium between cloud 
formation and destruction and has a large dynamic range of surface density, we find a fairly 
tight relationship between T^sfr and S, with 1 + p ~ 2 (see Fig. [T^ . This implies the 
fraction of dense gas follows M^^nsel Mtot = ^densel^ oc S. If we interpret this in terms of 
cloud formation/destruction equilibrium (cf. eq. [221) . with a constant mean cloud lifetime 
given by equation f|T2l) . this implies a dense gas formation time oc S~^. As discussed in 
§3.3.3[ our quantitative results are generally consistent with a formation time for dense gas 
(xtjoitff, which vary (exactly or approximately) oc S^^ in our models. 

In other recent numerical work, star formation prescriptions T^sfr oc ^/tjf have been 
adopted, where either all of the gas or just high-density gas is included in the right-hand 
side. This would imply T^sfr oc Tj^/'^{G / Hy/'^ for the dependence on surface density and 
disk thickness. For a disk in vertical hydrostatic equilibrium with vertical velocity dispersion 
the natural thickness varies as if oc cr^/(G'S), which would imply T^sfr oc T?G/az. Thus, 
a vertically-resolved disk with a constant vertical velocity dispersion would be expected to 
yield an index 1 + p = 2. If the disk thickness is determined not by hydrostatic equilibrium 
but in some other way, however, the resulting star formation rate and the index in the K-S 
law would depend on the numerical prescription (or physical process) that sets H. In our 
models, we have a flared disk H oc R and set S oc in our initial conditions, which 
accounts for the index 1 + p ~ 2 that we obtain. If, on the other hand, the value of H were 
constant in a given simulation (either by design for a two-dimensional simulation, or as a 
consequence of limited spatial resolution in a three-dimensional simulation), then the result 
would be 1 + p ~ 1.5. Thus, limited vertical resolution can potentially artificially reduce the 
scaling index in the K-S relation, as measured from numerical simulations. A fully-resolved 
vertical dimension is therefore required if the star formation prescription is to be based on 
a volume density. In practice, the resolution requirement can be quite demanding if the 
disk is dominated by cold atomic or molecular gas, since c^/(7rGE) = 4pc(T/100K)(S/10 
MqPc"^)"^. This also points to the necessity of incorporating turbulent processes in three- 
dimensional models, since observed cold gas is in fact dominated by turbulent rather than 
thermal pressure. If these turbulent effects were not included, the disk thickness would be 
unphysically small. 
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4.2. Model Limitations and Future Prospects 

4-2.1. Spiral Structure 

In spiral models, the external spiral potential is initially the primary driver for enhancing 
the density, leading eventually to the growth of clouds. In models that evolve for a significant 
amount of time, soon after feedback and the dispersal of cloud gas the global spiral pattern 
is disrupted, and eventually vanishes. With the simple feedback prescription that we have 
adopted, we were unable to simulate a spiral galaxy in which the global spiral pattern is 
maintained simultaneously as cloud gas is returned to the ISM through feedback. 

If the arms truly are long lasting, then either the spiral potential is much stronger 
than in our models {F » 10%), and/or the real feedback events are not as disruptive of 
structure on kpc scales . Very large F, how ever, does not appear consistent with observations 



of the old stellar disk (iRix fc Riekelll993l ). One possibility is that realistic feedback is both 
gentler and less correlated than the simple prescription of our current models, and as a 
consequence the spiral arm coherence would not be destroyed by large-scale shells. Indeed, 
semi-analytic models suggest that photo-ionization may evaporate much of the mass in a 
typic al GMC before the pr essure-driven expansion of HII regions unbinds the whole cloud 



le.g- 



Krumholz et al.l 120061 ). Those models do not include supernovae, however, which are 
unavoidable if a GMC survives for more than one generation of OB stars. Still, supernovae 
that are less correlated in space and time than the extreme case we have considered would 
disperse cloud gas in smaller parcels. Less correlated energy inputs would produce shells 
with diameters less than the spiral arm thickness, and could more easily leave global spiral 
structure intact. By studying how the resulting spiral morphology varies with the correlation 
of feedback energy, it will be possible to place limits on how correlated star formation is in 
real galaxies. 



4.2.2. Multiphase ISM 

The models discussed in this paper use the simplest possible prescription for gas ther- 
modynamics, which is an isothermal equation of state. Our adopted sound speed of c<j = 7 
km s~^ corresponds to a temperature of T ~ 10^ K, characteristic of the warm phase of 
the ISM. We adopted this approach in order to investigate, in a controlled fashion, various 
separate effects that can contribute to the regulation of star formation. 

In parallel with our simplified ISM thermodynamics, our approach to modeling feedback 
from star formation is also reduced to the most basic elements. In our models, we follow 



-24- 



the expansion of clouds subsequent to correlated SN events. Of course, in a real SN event, 
thermal energy is injected into the ISM, and it is the expansion of a very hot and very diffuse 
bubble of gas that drives the formation of a dense shell around it. The cooling time in the 
high density shell is short, so at late stages the isothermal approximation is adequate. The 
cooling time of the hot interior of each individual bubble, and of the hot phase of the ISM that 
results from merging SN remnants, is much longer. However, the hot phase contains only a 
very small fraction of the total ISM mass. From the point of view of most of the mass in the 
ISM, the primary effect of SN is to inject momentum. By adopting an isothermal equation 
of state, and treating feedback as providing momentum inputs, this effect is captured in an 
approximate way. 

A significant limitation of our models is that we do not treat the cold (T~100 K) 
atomic component of the ISM explicitly. Because the level of turbule nce in the atomic 



component is comparable to the thermal velocity dispersion of warm gas (IHeiles fc Troland 



20031 ). the effective pressure in the cold medium may be comparable to the thermal pressure 
in the warm medium. The dynamics associated with "turbulent pressure" may, however, be 
quite different from those resulting from micro-physical thermal pressure. A very important 
direction for future work is to study directly how large-scale gravitational instabilities and 
spiral structure develop in multiphase, turbulent, cloudy gas. 

Another limitation of our models is that they are two-dimensional (although the disk 
flares with radius). This constrains feedback energy to be confined within the galaxy's 
midplane, and does not allow for dynamically evolving disk thickness. In the real ISM, 
correlated SN may be important in driving the SN heated gas away fro m the midplane of 



the g alaxy into the halo, through so-called chimneys and superbubbles ( iNorman fc Ikeuchi 



19891 ). To explore the effect of this energy loss in an approximate way, in our models we 
consider both a "standard" momentum input per SN, and a momentum input reduced by a 
factor of ten. However, the cycling of gas through the galactic halo has other consequences 
as well. Afte r this gas is cooled in th e halo, it falls back onto the disk in the form of 



cloudlets (e.g. iJoung &: Mac Lowll2006l ). Even though r ecent simulations have shown tha t 



the fraction of mass that is vertically driven is small (e.g. Ide Avillez &: Breitschwerdtll2004l ). 
the in-falling clouds may still affect the dynamics of the disk and may also act as another 
source of turbulence. 

In order to accurately model disks that account for the effects of SN heating, chimneys, 
superbubbles, and the return of halo gas onto the disk, a three dimensional grid, as well 
as explicit treatment of heating and cooling, are necessary. Three dimensional simulations 
will also allow us to test the sensitivity of the Kennicutt-Schmidt slope to the disk thickness 
(which evolves in response to star formation), as discussed in § 14. 1[ These directions are 
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important avenues for future research. 

4.3. Summary 

In this paper, we consider the formation of self-gravitating structures in global models of 
spiral galaxies, focusing on the effects of star formation feedback. Our numerical simulations 
adopt a simple, isothermal treatment of the gas, and follow the flow in the disk by integrating 
the hydrodynamical equations on a polar grid. We incorporate vertical disk thickness effects 
within the solution of the Poisson equation, which assumes that the disk flares as H (x R. 
The feedback model treats the specific star formation rate in gas above a given surface 
density threshold as a constant, Rgj^fMsN- Feedback is implemented by spatially-resolved 
radial momentum injection subsequent to star formation events; the momentum injection 
is proportional to the number of stars formed. In order to explore the sensitivity of the 
resulting model properties to the feedback parameters, we consider a range of specific star 
formation rates, star formation efficiencies esF, and momentum injection per massive star 
Prad- We analyze the ISM spatial distribution, star formation rates, and turbulent properties 
of our model disks in cases with and without an externally-imposed spiral gravitational 
perturbation. 

Our main findings are as follows: 

(1) In models where Prad is comparable to the level expected from a supernova, clouds 
are destroyed by star formation events and the mean turbulence level is high. In models 
where Prad is a factor of ten lower, to represent inefficient feedback (e.g. if SN energy is 
vented vertically rather than kept in the disk), the self-gravitating structures that form 
are not destroyed by feedback, and the turbulence levels are substantially lower. Turbulence 
levels are insensitive to the star formation rate parameter Rsn and the overall star formation 
rate, however. 

(2) In models with strong feedback, expanding fiows lead to collisions of shells, which 
then lead to gravitational collapse of overdense regions and further star formation events. 
In this sense, our models are a concrete realization of the concept of self-propagating star 
formation. We find, however, that the net effect of feedback is to lower the rate of star 
formation. That is, when we compare models with strong feedback (large Prad) and weak 
feedback (small Prad), the former have lower resultant star formation rates. Similarly, when 
we compare models (at large and fixed Prad) that have high or low feedback event rates Rsn, 
the fraction of dense gas is lower when the event rate is higher. In principle, turbulence can 
either enhance collapse and star formation (by inducing shell collisions) or suppress collapse 
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and star formation (by breaking up overdense regions) . Our results show that although both 
effects occur, the latter dominates: star formation is in net suppressed by feedback. 

(3) For E ~ 10 — 100 M0 pc~^, the range in I^sfr for our simulations is similar to the 
range observed in normal disks. The slope of the Kennicutt- Schmidt scaling relation T^sfr 
E^"*"^ is steeper (1+p ~ 2) in our simulations than the slopes found from current observations 
at high (average) surface density. The discrepancy may be due to our assumption that the 
disk thickness varies with radius as H cc R. Indeed, our numerical results are consistent 
with the theoretical prediction that tsF tj or tff when the gravitational times tj and 
tff are calculated based on our model prescription. We point out that shallower scalings of 
^SFR with E would be expected if the vertical velocity dispersion increases with S. This 
would increase the disk thickness at small radii (where E is large) relative to what we have 
assumed, and consequently increase the gravitational times and reduce T^sfr- 

(4) Motivated by our own results, we remark that in general, the thickness of the gaseous 
disk in a galaxy (either observed or simulated) is important for setting the index in the 
Kennicutt-Schmidt relationship. Numerical simulations must resolve the natural disk scale 
height (set by pressure and turbulence) if the adopted prescription for star formation depends 
on the volume density p of gas. A simulation that is vertically unresolved {H const.) 
while adopting psF oc p/tff{p), and hence Hsf oc ^/^//(p)) automatically yield an 
index 1 + p = 1.5 in the K-S law since tjj oc (T./H)^'^. Fundamental understanding of K-S 
laws requires a self-consistent determination of the dependence of on E. 

(5) For turbulence driven by expanding shells in overdense regions, we find that the 
power spectra decrease with decreasing size consistent with the scalings for shock-dominated 
flows ("Burgers turbulence"). While typical mass- weighted velocity dispersions on kpc scales 
in our high-Prad models are 31 km s^^, these decrease to 10 km s^^ on 100 pc scales. Radial 
and azimuthal components of the velocity dispersion in a given scale are comparable. 

(6) For all of our models, the maximum masses of dense clouds that form are several 
million Mq, consistent with observations of the upper cutoff in GMC/GMA mass distribu- 
tions in local group galaxies. In models with strong turbulence, such that self-gravitating 
condensations can form in colliding flows, a wider range of cloud masses results, with a lower 
peak in the distribution (but similar upper cutoff). Higher resolution simulations will allow 
for a more detailed analysis of the mass distributions. 

(7) Within the context of the feedback prescription and parameters for our current set 
of models, we flnd that global spiral patterns are not maintained. For low Prad, insufficient 
momentum is injected to overdense structures so that arm clouds continue to collapse, even- 
tually depleting the surrounding spiral arm gas. For high Prad, large-scale expanding shells 
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form and the global spiral structure is destroyed as cloud gas is dispersed. We conclude 
that highly-correlated star formation, which is the limit that we adopt in the present mod- 
els, is incompatible with long-lived spiral structure. It will be interesting to determine, by 
comparing spiral morphology with results from models adopting differing feedback prescrip- 
tions, what constraints are placed on the spatial and temporal correlation of star formation 
feedback in real galaxies. 
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A. Appendix 



In the Appendix of Paper I, we described two methods to solve Poisson's equation 
numerically on a polar grid; both methods employ Fast Fourier Transforms (F FTs). One 
method sums the potential from concentric rings, as described by iMillerl ( 1l976l ). The other 
method employs a coordinate transformation from polar coordinates to a Cartesian-like 
coordinate system. The former method is exact, but computationally expensive, and the 
latter is an approximation, but computationally efficient. 

H ere, w e describe another FFT based method that is exact , and more efficient than the 

Milleil fll976h methodQ The basic scheme is described in lKalnajsl ( 119711 ) and lBinney fc Tremaind 
(119871 ): we describe a modification of Kalnajs's method that includes the effect of nonzero 
disk thickness H, which also acts as softening. 



The potential $ at each position {R, 0) on the disk, at z=0, is 



J J J [R'^ + - 2RR' cos{(P - (j)') + z'^] 2 

Here, G is the usual gravitational constant, S is the total surface density, and the function 
/ = p{z', R', (j)')/H{R', (j)') describing the vertical profile of the volume density must be nor- 
malized, J^^dz' f{z' , R' , (J)') = 1. Substituting u' = Ini?', and (' = z' / ^/2R' in equation 



^As in Paper I, we again make use of the freely-available FFTW software (jFrigo fc Johnson 



20051) . 
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f lAll) . the potential reduces to 

HR, 0, . ^ 0) ^ -Ge^ [ du' [ d^' [ dC e--^e-nCW,mi<^') ^ 

J J J [e«-«'(cosh(ti-«') -cos(0-0')) + C'']^ 

If R'p{z',R',(j)')/E{R',(f)') = e-^' f{C,u',(j)') = g{C) is a function of C' only (see below), we 
can define 

I(u' -u,(P'-<p) = e"'-" [ dC ^ r- (A3) 



Using the definition of / in equation (lA2p . we obtain $ as a two dimensional convolution: 
^{R, ^,z = 0) = -Ge" j du' j d(P'T.{u\ 4>')I{u' -u,<p' -<p). (A4) 



Applying the Fourier convolution theorem to equation flA4l) . the gravitational potential can 
be computed by taking the Fourier transform of S to obtain S, and then taking the inverse 
Fourier transform of the product of S and /, where J is the Fourier transform of /. In 
hydrodynamic simulations, / can be computed once at the beginning of the simulation run, 
so that only two FFTs need to be performed at each timestep, FFT(S) and FFT^^(S/). 

The function J, and therefore its convolution J, depends on the normalized vertical 
distribution function g{C)- For the specific case of a Gaussian vertical density distribution 
(which holds if the vertical gravity is dominated by that of the stellar disk), 

f{z,R) = — (A5) 

For a disk that flares as H{R) oc R, we define H = H{R)/R, so that 

e"/(C,i?) = ^=- = ^7(C)- (A6) 



Similarly, if the vertical density follows a sech^ distribution (true if the gaseous self-gravity 
dominates), then ^(C) = i2n)-hech'^{CV2/H). 

For our simulations, we adopt the Gaussian profile; this yields the following explicit 
expression for /: 

I(u' -u,<f)' ~<f)) = ^^ dC -. (A7) 

V2^n J-oo [e"-«'(cosh(M - u') - cos(0 - 0')) + C'']5 



Finally, we comment on the assumption H[R) oc R which enables the three-dimensional 
gravitational integral to be written as a two-dimensional convolution. If the stellar disk 
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dominates gravity, then for an isothermal disk the vertical density distribution is Gaussian 
with H/R = CsQ^{c^^z/c*^r)/2vc, so values of Cg/vc, Q*, and c^^^/c^^/j that are independent of 
radius imply constant H/R. Similarly, if gas is the dominant component for vertical gravity, 
H = cI/tcGT,, so that H/R = CsQ/\/2vc. If both the Toomre Q parameter and Cg/vc are 
independent of R, then H/R = constant. For self-gravitating gaseous disks, if Q = 1 and 
t^r/c o=30, then H/R = 0.02. Including stellar gravity typically reduces if by a factor of ~ 2 



(e.g 



Kim et al.ll2002l ). 



For the simulations described in this paper, we use H/R = Ti = 0.01 in equation 
(IA7I) . We have tested other values of 7i, and find that our results are not sensitive to the 
exact value. However, large changes significantly affect the rate of growth of self-gravitating 
perturbations. 
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Pig. 1. — Density snapshots of Qo = 1 non-spiral model before any feedback, at time t/tcn-b 
— 0.84. Gray scale is in units of log(E/Eo). 
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Table 1. Parameters for Standard^ Models 



Model 


Qo 


m 


RsN° 




Prad (10' Mq km s-i) 


(1) 


(2) 


(3) 


(4) 


(5) 


(6) 


QIA 


1 





1 


0.025 


0.34 


QlB 


1 





1 


0.05 


0.34 


QlD 


1 





10 


0.05 


3.4 


QISA 


1 


4 


1 


0.025 


0.34 


QISB 


1 


4 


1 


0.05 


0.34 


QISC 


1 


4 


1 


0.05 


3.4 


QISD 


1 


4 


10 


0.05 


3.4 


QISE 


1 


4 


10 


0.05 


0.34 


Q2SA 


2 


4 


1 


0.025 


0.34 


Q2SB 


2 


4 


1 


0.05 


0.34 


Q2SC 


2 


4 


1 


0.05 


3.4 


Q2SD 


2 


4 


10 


0.05 


3.4 



^1024x1024 zones; R e 4-11 kpc; (j) e 0-f radians 

^Units of 2 X 10~-^^Mq^ yr~^, i.e. number of SN per 50 years per 
lO^Mn 



Fig. 2. — Density snapshots of model QlA, at time t/tori — 1.125. Gray scale is in units of 
log(E/Eo). 




Fig. 4. — Qq = 1 spiral model, without feedback, at t/torh = 0.675. Gray scale is in units of 
log(S/So). 
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Fig. 5. — Qq = 2 spiral model, without feedback, at t/torh = 1-04. Gray scale is in units of 
log(S/Eo). 
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Fig. 6. — Models QISC (left) and QISD (right) at t/torb — 0.73. Gray scales are in units of 
log(E/Eo). 
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Fig. 8. — Cloud masses in models with strong (thick lines) and weak (thin lines) feedback. 
Left: Histogram of Md in models QID (thick) and QlA (thin), up until time t/torb = 1-125 
(see Figs. [2]-|3]). The mean (median) M^i for models QlA and QID are 1.2x10^ (0.8 xlO^) 
and 0.7x10^ (0.5x10^) Mq, respectively. Right: Histo gram of M^i in models QISD (thick) 
and QISA (thin), up until time t/torb = 0.73 (model QISD is shown in Fig. [H]). The mean 
(median) M^i for models QISA and QISD are 2.2x10^ (1.9x10^) and 0.8x10^ (0.6x10*^) 
Mq, respectively. 
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Fig. 9. — SFR (top) and Vturh (bottom) for models with (5o=l; without spiral structure. 
The values in parentheses in the legend are the SN rate parameter Rsni the star formation 
efficiency esF, and SN momentum Prad (in M© km s~^) of each model. The large open squares 
in the bottom panel are the turbulent velocities for a simulation without any feedback. 
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Fig. 10. — SFR (top) and Vturb (bottom) for models with Qo=i, as in Figure [9l but with 
spiral structure. 
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Fig. 11. — SFR (top) and Vturb (bottom), as in Figure fTOl but for models with Qo=2, with 
spiral structure. 
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Fig. 12. — Mass weighted turbulent velocities vs. mean surface density of QIS models, 
averaged in annuli of widths 1 kpc, and in the time interval t/torb £ 100 — 116 Myr. Turbulent 
velocities are only shown from annuli and time intervals within which feedback events have 
occurred. 
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Fig. 13. — Turbulent power spectra of Model QID. Power is shown at constant kji (left) and 
constant A;,^ (right) (each slice is along the minimum nonzero value of the respective k). To 
obtain the dimensions of k(p {k^p = mn^p/R), we use the mean radius of the grid. Best fits for 
values between log(A;) e — 1 gives slopes of -2.9 (left) and -2.4 (right). 
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Fig. 14. — Schmidt law for models in Table [H Each point is obtained by binning the 
simulation data in radii, with 1 kpc widths, and in time, with 18 Myr widths. Only annuli 
and times with at least 1 feedback event are included. The lines are the best fit to the points; 
with their slopes ranging from 1.8 - 3.8 (top), 1.1 - 3.0 (middle), and -2.0 - 1.9 (bottom). 
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Fig. 15. — Left: Star formation rate vs. surface density for model QID (triangles; i? e 4 - 
11 kpc), as well as the corresponding model of the inner region (crosses; R G 0.8 - 2.2 kpc). 
Best fit lines for each model are also shown, with slopes of 2.4 for the 4-11 kpc model 
and 2.2 for the 0.8 - 2.2 kpc model. Right: Triangles and cros ses from fi g ure on the left 
are shown, along with globally averaged observational data from iKennicuttI ( jl998bl ): circles 
show normal spirals, with best fit slope of 1.3, and diamonds show IR starburst sources, with 
best fit slope of 1.4. 
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Fig. 16. — Star formation times from model QID, as a function of mean surface density. 
Different symbols correspond to different time bins, of width 18 Myr. Simulation data are 
also binned in radii with widths 1 kpc. Filled symbols show actual depletion times (eq. |20j). 
and open symbols show predicted depletion times (eq. [21]), for each annular and temporal 
bin. 
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Fig. 17. — Star formation times from model QID, as a function of Jeans time (top) and 
free-fall time (bottom). Simulation data are binned as described in the caption to Figure [T6| 
with the innermost and outermost annuli excluded. The dashed lines, shown for comparison, 
have slopes of 1. 



